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. ! ABSTRACT 

o . 

•7- , We explore the possibility that massive black holes comprise a significant fraction 

of the dark matter of our galaxy by studying the dissolution of galactic globular clusters 
bombarded by them. In our simulations, we evolve the clusters along a sequence of 
King models determined by changes of state resulting from collisions with the black 
holes. We include mass loss in collisions as well as the heating of the remaining bound 
■ stars, and determine the role which a finite number of stars plays in the variance of 

' the energy input and mass loss. Several methods are used to determine the range of 

black hole masses and abundances excluded by survival of galactic globular clusters: 
simple order of magnitide estimates; collision-by-collision simulations of the energy 
input and mass loss of a stellar cluster; and a 'smoothed' Monte Carlo calculation of 
] the evolution of cluster energy and mass. The results divide naturally into regimes 

ON ' of 'small' and 'large' black hole mass. 'Small' black holes do not destroy clusters in 

(-h \ single collisions; their effect is primarily cumulative, leading to a relation between Mbh 

O and /halo, the fraction of the halo in black holes of mass Mbh, which is /haioMbh < 

constant (up to logarithmic corrections). For /halo = 1 , we find Mbh ^ 10 3 M Q by 
requiring survival of the same clusters studied by Moore ( 1993| ), who neglected cluster 



o 

r/^ \ evolution, mass loss, and stochasticity of energy inputs in his estimates, but reached a 

similar conclusion. 'Large' black holes may not penetrate a cluster without disrupting 
it; their effect is mainly catastrophic (close collisions), but also partly cumulative 
(distant collisions). In the large Af D h limit, /h a i (but not M D h) can be constrained by 
computing the probability that a cluster survives a combination of close, destructive 
■ encounters and distant, nondestructive encounters. We find that it is unlikely that 

/halo 0.3 by requiring 50 per cent survival probability for Moore's clusters over 10 10 
years. 

Key words: Galaxy: general - globular clusters: general - Galaxy: halo - Galaxy: 
kinematics and dynamics - dark matter 



1 INTRODUCTION 



In this paper, we reexamine the idea that black holes compose a substantial fraction of the dark matter in the halo of our 



galaxy. Ostriker & Lacey ( 1985) originally proposed that heating by a significant population of halo black holes with masses 

A^bh ~ 10 6 Mq could explain Wielen's (1977) inference that the velocity dispersions <t* of disk stars of age t* follow the trend 

1/2 

cr* oc tj . Although subsequent work questioned the basis of this argument, both because it is unclear that cr* rises as fast as 



1/2 



(e.g. Carlberg et al. 1985; Stromgren 1987; Gomez et al. 1990), and that black holes are needed to explain the observed 
trend of <r* with t+ (e.g. Lacey 1991), the heating of disk stars remains a constraint on the mas ses and abundance of any 
hypothetical population of massive halo objects; in the notation of Wasserm an fc Salpeter ( 1994 ), if the halo consists of a 
fraction /halo in the form of massive black holes, then the Ostriker & Lacey (1981:) argument implies /haioMbh ^ 1O 6 ' 6 M0. 



Limits on the masses and abundances of massive black holes could have important cosmological implications (e.g. Loeb 1993; 
Umemura, Loeb, & Turner, 1993; Loeb & Rasio 1994). 
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Carr (1994) has reviewed various limits on baryonic or black hole dark matter in the halo of our Galaxy. (See also Carr 
& Sakellariadou 1998.) One of the most powerful constraints on black hole properties yet proposed was put forth by Moore 
(1993) who concluded that the survival of a set of relatively tenuous, low mass (M < 4 x 10 4 Mq) globular clusters over a 
timespan >7 x 10 9 years would be possible only if Mbh ^ 1O 3 M0. The argument employed by Moore (1993) was first applied 
to this problem by Wielen ( 1985| ), who explored perturbations of more massive clusters (M ~ 10 6 Mq) by heavier black holes, 
in the range advocated by Ostriker & Lacey (1985). [Klesson & Burkert (1995) extended Moore's and Wielen's calculations 
to a range of globular cluster properties.] The idea is to compute the heating of a globular cluster by passing black holes 
over a chosen timescale comparable to the cluster age [as noted above, Moore (1993) used 7 x 10 9 years, roughly half the age 
inferred for typical clusters]. If the total energy imparted by black hole perturbations is sufficiently large, then the cluster is 
said to be disrupted. 

In view of the importance of this problem, we have begun a more comprehensive study of the disruption of the same 
set of globular clusters considered by Moore (1993) by a hypothetical population of halo black holes. This investigation aims 
to tighten up Moore's argument in several different ways. Some of the improvements are technical [e.g. Moore used a simple 
analytic form for the energy input to a cluster by a passing black hole which is only very approximate; Klesson & Burkert 
(1995) improved on his formula], but others are qualitative. Of particular importance are: (1) Moore computed the energy 
input for a 'static' cluster, whose structure was held fixed. Our calculations evolve the clusters along a King sequence. (2) 
Qualitatively, one expects a cluster that gains a large amount of energy compared to its initial binding energy from encounters 
with passing black holes to lose a large fraction of its mass, too. Our calculations include mass loss by the clusters. (3) Although 
one can get a rough idea of the survival probability by considering the mean heating of a cluster, the energy transfer process 
is actually stochastic, and the variance in the energy input may play an important role in final estimates of critical masses for 
disruption to be likely. (4) Moore's calculations pertain to black holes of relatively low mass, which are incapable of disrupting 
a cluster in a single perturbation. For black hole masses Mbh MV Te i/a c i, where <t c i is the characteristic velocity dispersion 
of the cluster and V Te \ its characteristic speed relative to the approaching black hole (V Te \ ~ 330km s _1 is a fair estimate), 
a single encounter at the cluster's tidal radius will likely destroy it. In this regime, one can obtain limits on the halo mass 
density in massive black holes, but not on Mbh (e.g. Wielen 1988; Wasserman & Salpeter 1994). 

Our study combines analytic and Monte Carlo calculations. Our most complete results come from Monte Carlo simulations 
in which we simulate each encounter between a cluster and a black hole separately. For simplicity, we shall consider clusters 
at fixed galactocentric radius, as was done by Moore (1993). We model clusters using N point masses whose positions and 
velocities are chosen from a King distribution (see, e.g., Binney & Tremaine 1987). Velocity perturbations are computed for 
a given black hole impact parameter and speed relative to the cluster center of mass using the impulse approximation (e.g. 
Binney & Tremaine 1987). The change in velocity of the cluster center of mass is subtracted from each individual velocity 
perturbation to determine the change in cluster energy and mass as a consequence of the collision. Determining which stars 
are ejected is tricky in any scheme that does not employ a direct N body simulation of interparticle interactions, but as long as 
the perturbations in individual encounters are not excessive, it should be sufficient to designate for ejection those stars whose 
post-collision velocities exceed their local pre-collision escape speed. To find the new King model that describes the remaining 
cluster, we need its tidal radius in addition to its mass and energy after the black hole encounter; we get this by assuming that 
the tidal radius is proportional to M 1 ^ 3 , consistent with our assumption of fixed galactocentric radius. Within these 'rules 
of the game' we simulate the evolution of the globular clusters studied by Moore (1993) over a timespan of 10 10 years, or 
until they are disrupted, whichever comes first. Our simulations are more comprehensive than the Monte Carlo calculations of 
Klesson & Burkert (1995), who chose discrete black hole encounter times, relative velocities V fe \ and impact parameters b from 
the appropriate probability distributions, but merely updated the cluster velocity dispersion by a completely deterministic 
amount Aa(b, V Te \), without accounting for stochasticity of heating, mass loss or the change in internal cluster structure in 
individual encounters. 

Several different criteria are employed to decide whether or not a cluster has been destroyed. Some of these may be called 
'global' in that they depend on integrated properties of the cluster, such as total energy or total mass. We can regard a cluster 
as having been destroyed, for example, when its mass or energy or energy per mass has changed by a fractional amount in 
excess of some pre-set values (e.g. 0.5). Since we evolve models along a sequence of King models which is limited, clusters may 
also die when they reach the end of the sequence (see also Chernoff, Kochanek, & Shapiro 1986). The other criteria for cluster 
disruption to be used are 'local', in that they depend on changes in the properties of a cluster as a consequence of a single 
collision. Thus, if the mass or energy of a cluster changes by more than some pre-set fractional amounts in an encounter with 
a black hole, we shall regard it as disrupted. This should also allow us to control the inaccuracy of our criteria for determining 
the mass loss per encounter somewhat. We evaluate survival probabilities for the various criteria separately. 

The Monte Carlo calculations outlined above allow us to study both the large and small Mbh regimes. In the large Mbh 
regime, where destruction occurs only after relatively few encounters, the Monte Carlo calculations ought to be reasonably fast 
computationally. However, for small Mbh where the encounters are more frequent (the regime focussed on by Moore 1993), 
we expect the Monte Carlo calculations to be more cumbersome, thus limiting the number of stars we can use in realizing 
clusters. Ideally, one would like to be able to represent the stars 'one-by-one' since the variances in energy input depend on 
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the number of cluster particles. Fortunately, the perturbations due to individual encounters are relatively gentle in the small 
Mbh regime, and the problem lends itself to a Fokker-Planck treatment, which promises to be faster (at the price - justified 
in our view - of losing the ability to resolve short timescale structure in the cluster evolution). We present a two dimensional 
Fokker-Planck scheme in which we follow cluster evolution in mass and energy. 

The calculations reported here suffer from two or three principal deficiencies. Most important is their reliance on the 
sequence of King models and on the impulse approximation. In addition, one would like to interweave perturbations by a 
hypothetical population of black holes with well-established sources of heating, such as disk shocking (e.g. Ostriker, Spitzer, 
& Chevalier 1972; Spitzer & Chevalier 1973; Chernoff et al. 1986; Binney & Tremaine 1987); in the similar problem of wide 
binary evolution, the interplay of perturbations by stars, molecular clouds, and dark matter is known to be important (e.g. 
Retterer & King 1982; Bahcall, Hut, & Tremaine 1985; Weinberg, Shapiro, & Wasserman 1987; Wasserman & Weinberg 1991). 
Moroever, globular clusters evolve on their own as a consequence of internal relaxation, resulting in evaporation and energy 
changes even without external perturbations; limits on halo black hole properties may be altered when internally induced 
changes in state are accounted for properly. These important effects will be ignored here to concentrate solely on the cluster 
evolution due to the collisions with black holes. We shall study some of these issues in a subsequent paper (Murali et al. 1998). 



2 QUALITATIVE OVERVIEW OF THE COLLISION PROCESS AND EVOLUTIONARY SCENARIOS 
2.1 Approximations 

Throughout this paper we employ the impulse approximation and ignore deflection of the perturbing black hole orbit from a 
straight line. For nonpenetrating encounters at impact parameter b and relative velocity V ro \, these approximations are valid 
when D.b/V vc i < 1. The characteristic frequency is roughly Q, ~ (GM/r^) 1 ^ 2 , where M is the cluster mass and 9"h m is the 
cluster half mass radius, for the outer parts of the cluster where the tidal perturbation is strongest. If we set V Ye i = £i V c , where 
V c is the galactic circular speed, and let ao = ^(GM/rhm) 1 ^ 2 be the cluster's central one-dimensional velocity dispersion, 
then the impulse and straight line approximations should hold for 6/ri lm < (£,i£,2)V c /<to- For the clusters studied in this paper, 
Vc/(Tq ^ 10 — 100, and the approximations only fail far outside ri lm . 

For penetrating encounters, this assessment remains valid for 'typical' collisions, but may fail for perturbations of particles 
deep in the cluster core. This is because the characteristic frequency near the center of the cluster where the perturbations are 
largest is fl ~ (47rGpo/3) 1,/2 , where po is the central mass density of the cluster. Adopting b = r^m as typical, we find that the 
impulse and straight line approximations hold for Kci ^ 0.1 x r hm (pc)[po/M0pc _3 ] 1 ' /2 km s _1 . Although our approximations 
would fail for the most concentrated clusters, they remain true for the relatively tenous ones studied here (and indeed for 
many 'normal' clusters). 

We also ignore all processes influencing cluster evolution except perturbations by black holes even though our limits are 
based on survival probabilities over timescales long enough for tidal shocking and internal relaxation to be important. In more 
realistic simulations, including these effects could tighten limits on properties of hypothetical halo black holes. 



2.2 Single Collisions 

Consider a single collision between a globular cluster and a black hole which takes place in the halo of the Galaxy. Let the 
cluster have tidal radius r t , N stars, King model parameter Wo = ipo, and total energy E. The black hole passes the cluster 
with a speed V Te \ and an impact parameter b. Define the 'collision parameter', rj c — (Mbh/M)(cro/V r rc i). In section §^^, it will 
be shown that in the impulsive limit the energy input, mass loss and their variances are 



AE(b,V le l) _ ■ 



\E\ 



rj c R E {b, rj c , ipo) 



AM(b, V Icl ) a - 
= -r) c RM{b,rjc, i>o) 

-ga" = jjREE(b,rj c ,ipo) 

2 2 

-^f" = JfRMAt(b, Tj c , Ipo) (1) 

where the b dependent radial functions have the approximate limits (see §p.3|) 



R E ~ C E ^ ) 
Rm^Cm^o) 
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Ree^C ee (4> ) 

Rmm ^ C M m$o) (jf) (2) 

when b > r t . As much of the mass in the cluster is contained within the core radius, one can actually use the b~ 4 dependence 
all the way into the core as a first approximation. Examples of the radial functions Re, Rm, Ree, and Rmm for impacts 



inside the cluster are given in §5.2 



For large enough Mbh, a single collision will suffice to disrupt the cluster. This will be referred to as the 'high black hole 
mass limit'. Define 'disruption' to occur for an energy input of size f\E\. Then the cluster is destroyed in a single collision for 

- < W^Y'V. (3) 

n n \ f J 

A safe overestimate of the mass, Mhi g h, at which the cluster may be disrupted is given by the black hole mass at which bd = r t . 
Using V [c i = £iV"c, £i — 1.5, and the Galactic circular speed V c = 220km s _1 , the result is 

m »+ = m vMT (4) 

which is bigger than M by a factor of V c /o~o S> 1. 

The probability distribution for energy input (and mass loss) can be understood in the following qualitative terms. The 
ratio of the variance to the mean energy input is 



°"ab Ree Cee 1 b 4 



(AE) 2 NrgB* NC%rf c r$ 
This becomes small for large N and small Mbh, but is also proportional to b 4 , insuring that ctab > Ai5 for 



(5) 



a 



b_ > bdiff _ / 
n ~ n V Cee 



1/4 



N 1/4 V 1/2 . (6) 



Note that b^is > r t for Mbh/M > V le \/(VNao)- The distributions of mass loss and energy input are sharply defined about the 
mean for b < fediff; and become 'fuzzy' for b > &diff- In particular, the energy input will always be sharply peaked about the 
mean at the destructive radius fed, because &d/bdiff = (Cee / /CeN) 1 ^ 4 <C 1. Hence the probability of getting an energy input 
which is not destructive when b < bd is quite small. 



2.3 Evolution Over Many Collisions 

If Afbh <§C Mhigh ('the small Afbh limit'), it will take many collisions to disrupt the cluster. Individual collisions only 'tickle' 



the cluster, and the evolution can be approximated by averaging over the effects of many collisions. In §5.4, we find that the 
mean rate of energy input, averaged over 6 and V rc i, takes the form 

@ = d E (Vco,ipo)r (7) 

where F = n bh Trr^V c , n bh is the number density of black holes, rj c0 = (M hh /M)(ao/V c ) and d E {ri c0 , ip ) ~ k £ (^o) hi(l /r} c0 )rjl . 
The mean time for disruption, which we define here to occur at AE = f\E\, is 

T - _ / 1 rs s 

disrupt " (E) ~ r « B ln(l/^ )^ ^ M bh [ ) 

so that for very small Mbh it takes a long time to disrupt the cluster. Given a value for Tdisrupt there is a certain critical value 
of Afbh, called Mbh,crit> above which the cluster is disrupted. The variance of the averaged energy input in a time T takes on 



the form (see §5.4) 

a AiS _ d EE . x 

~~EP~ - — ToT (9) 

where dEE^d),^) ^ «:eb(^o) ln(l/jj c0 )?jcO- 

Given an ensemble of clusters with initial energy Eq, a time T S> 1/ro later the cluster energies will be roughly E(T) — 
Eq + (AE(T)) ± o~e{T). For small times the variance will dominate the mean corresponding to a very broad spread in cluster 
states but for sufficiently long times the ensemble will peak sharply about the mean. We find that for times T ~ Tdisrupt 

^AE^disrupt) KEE \E\ KEE 1 „ , , 1rl , 

oc — < 1 (10) 



((AS^disrupt))) 2 Nk e (AE(T disIupt )) NfKE N 
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Figure 1 . Schematic of fraction of cluster destroyed as a function of Mbh 



so that the range of final states is always sharply defined about the mean for times long enough to disrupt the cluster 
when Mbh <C Mhigh. The probability of survival, P s , is given (here) by the fraction of clusters with AE < f\Eo\. This 
fraction will change from unity to a minimum value in a 'transition region' with size, <5Mbh,crit, dictated by the ratio 
(TAJ5(2di 8 ru P t)/(AS(TdiBrupt)> SO that <5A/ b h,orit/M b h,orit ~ TV" 1 / 2 < 1. 

When Mbh ?J Mhigh, any collision inside the destructive radius fed will destroy the cluster. In the tidal limit, the expected 
number of destructive encounters, Nd, in the time T is 



v = , t) 1/2 m 7 "^ 



(ii) 



where pbh = Mbb^bh is the mass density in black holes. For an ensemble of clusters with initial energy Eq, the destructive 
b < bd collisions act as a 'sink' for clusters while the nondestructive 6 > fed encounters give rise to a slower, diffusive energy 
change. The probability of survival, P s , (which takes on its minimum value in the Mbh ^ Mhigh limit) is the product of 
exp(— Nd), the probability that no single destructive encounters occur, and the probability that all the b > fed collisions 
combined give AE < f\Eo\. Since Nd does not depend on Mbh, only pbh (or equivalently /halo) can be constrained in the 
Mbh ^ Mhigh limit, not the black hole mass. In addition, the cluster evolution is quite stochastic in this regime as it depends 
on whether the destructive collisions do or do not occur. 

The model for cluster evolution described in this section is summarized in the (purely illustrative) diagram of Fig. |l[ 
The fraction of clusters destroyed in time T, called /dest = 1 — P«, changes from to ~ 0.8 in a region of width ~ 2OOOM0 
centered on Mbh.crit ~ 5OOOM0. Note that /d es t in this example does not asymptote at one, but instead reaches /dost ~ 0.8 
corresponding to Nd ~ 1. 
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Table 1. Parameters of the Globular Clusters 



Cluster Name M rt r COIC R g 

(M e ) (pc) (pc) (kpc) 



AM 4 700 11.7 3.7 30.0 

Arp 2 18000 75.0 9.5 20.4 

NGC 5053 37700 72.0 10.9 16.7 

NGC 7492 10400 45.4 4.5 18.7 

Pal 4 24900 92.3 15.4 96.0 

Pal 5 13700 76.4 13.9 16.5 

Pal 13 3000 30.0 3.0 25.7 

Pal 14 10400 117.5 22.4 69.9 

Pal 15 15000 41.9 10.5 30.0 



3 THE CLUSTER MODEL AND EVOLUTION OF THE CLUSTERS 

The nondimensional King models (e.g. Binney & Tremaine 1987) are uniquely determined by the normalized central potential 



ipo = = O)/^! = 9o/af, which is the parameter Wo of King (196C). The distribution function is given by 



A ing (r lW )d 3 rd 3 * = dW.^ i7¥ (e^^ 2 )/^)-l). (12) 

for stellar positions r, velocities v, and local escape speed v CS c(r). The dimensional King models can be specified by three 
independent quantities such as E, M, and r t . Define i>(V>o) by E = v(ipo) (GM 2 /n). For 0.0 < ipo < 8.5 , or —0.60 > v > 
—2.13, v(tpo) is single- valued so that the King model is known given E, M, and r t . A physical reason for excluding large ipo is 



that simulations have shown King models become susceptible to gravothermal instability at ipo ^ 7.40 (Wiyanto et al. 1985). 
The clusters we study in this paper are not core collapsed. 

A further restriction on the clusters is that they be tidally limited by the galaxy. To include the time-dependent effect of 
tidal stripping due to the galaxy would require detailed restricted three-body simulations for a range of both globular cluster 
orbits and orbits of stars in the clusters. In this paper we consider a first approximation in which the the galactic tidal field 
provides a relationship between rt and M, but does not contribute a time-dependent perturbing force. For circular orbits, if 
we define n as the distance from the cluster center to the Lagrange point of the cluster plus galaxy potential, we get 

i£ ~ Mi (13) 

for a point mass galaxy with mass M g and galactocentric radius R g , and we adopt M/rf = constant for our clusters as they 
evolve due to collisions with black holes. 

A weakness of our paper is its dependence on the King model sequence. As we shall see in for some clusters black 
hole collisions may force ipo — > after only modest energy input and mass loss, leading to very tight (but somewhat artificial) 
bounds on black hole properties. We shall rectify this deficiency in a subsequent paper where cluster structure is not restricted 



to the King sequence (Murali et al. 1995) 



We shall use two different sets of globular clusters to determine a maximum allowed black hole mass, Mbh,crit- First we 



examine the set of loosely bound globular clusters found in Moore (1993) and listed in Table 1, but then we also investigate 



a larger cluster perhaps more representative of the initial cluster population. 



4 THE BLACK HOLE MODEL 

Black holes of mass Mbh are assumed to compose a spherical halo with an isotropic velocity distribution. A fraction /halo of 
the total halo mass is presumed to be in the black holes. When no value of /halo is explicitly stated, /halo = 1 is assumed. 
We model the black hole population as a singular isothermal sphere with one-dimensional velocity dispersion Obh and mass 
density 

Pbh(Rg) = /halo 2 % QJ{2 

(14) 

(e.g. Binney & Tremaine 1987, Ch.4) with <Tbh = Vc/%/2 = 220/\/2km s _1 ~ 156km s _1 . The number density of black holes 
at R g is 

V c 2 

n hh (R g ) = /halo . — 52 ■ (I 5 ) 
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We ignore rotation of the black hole halo, so their velocity distribution is /(Vbh) = ( 7I V C 2 )~ 3/ ' 2 exp (— V 2 h /V 2 ^ . For 
computing the impulsive mass loss and energy input to the cluster due to the collision, we need the distribution of relative 
speeds, V re i. Let V c ] be the cluster velocity, so V ro i = Vbh — V c i. The distribution of relative speeds between the halo of black 
holes and the cluster becomes 

f(V lel )dV iei = ^^(exp^ - 2 j-exp(^ - 2 J)-^- (16) 

In this paper the clusters are on circular orbits with V c \ — V c . 
The mean relative speed given by this distribution is 

(Kd) = V j(z + ^)erf(:r) + -i=exp(-z 2 )J (17) 

where x = Vi/V and erf(s) = J* exp(— t 2 )dt. As x — > 0, {V re i) — > 2V c /\/tt (the Gaussian result), as x — > oo, 

(Vrel) —> Vci, and for Ki = V we find {V TC i(x = 1))/V e = £1 ~ 1.47. Note that the distribution in eq. ( |l6| ) is different from that 
used by Moore ( [l993| ) and Klesson & Burkert fll995| ). In addition, they chose to approximate the differential rate of collisions 
dT = nhh2nbdbf(V rc i)VrcidV ra i by nbh27r&db/( Vci) VdVci, which will lead to errors in the number of collisions and the rate of 
energy input and mass loss to the cluster. 

For a globular cluster at R g , there will be a certain value of Mbh below which there will be more than one black hole 
inside the cluster on average at any time (the 'many-body' limit). The number of black holes inside a cluster is A^ na i c i c = 
n bh 47rr t 3 /3 = (V c 2 r'i)/(3GM bh R 2 g ). Since r t ~ GM/al eq. [jj and Q give 3(V 2 / R 2 g )(r 2 /a 2 ) ~ 1 and so N insidc ~ M/(10M bh ). 
This exceeds one for Mbh < 0.1M which we will see is close to black hole mass limits for some clusters. The mathematical 
description of the energy input is complicated in this regime since the duration of a collision ~ rhm/Kci is longer than the 
time between collisions ~ (nbh7rrb m Kci) _1 . 



5 THE IMPULSIVE ENERGY INPUT FOR A SINGLE ENCOUNTER 



A» = — -— (18) 



In the impulse and straight line approximations, the velocity change of a star due to the passage of a black hole is 
2GM bh s 

VrclS 2 

where s _L V ro i is the projected vector from the black hole to the star at closest approach. Since the stellar velocity « < V c i 
and v <C V b h, we may neglect v in the relative velocity so that V re i — Vbh — V c i- In this approximation, all stars recieve a 
velocity kick in the same plane perpendicular to V fe i- The density of cluster stars can be projected onto this plane; then for 
a star at projected position R relative to the cluster center, s = R b, where b is the impact parameter of the black hole 
relative to the cluster center. 

The cluster is destroyed for impacts with Av > ao ■ An order of magnitude estimate of the ratio of these two speeds is 

A?; Mbh o-Q rf _ , , 

a ~ M Vd & 4 ~ nc & ( > 

For a penetrating impact with b < rt, the cluster will be destroyed if r\ c > 1. As the velocity kick, energy input, etc. in the 
impulsive limit must scale as Mbh/Vci, r\c is a convenient dimensionless measure of the destructiveness of the collision. 

Next let star i have a mass m; and a velocity Vi with respect to the center of mass before the collision. The energy of 
the cluster before the collision is then 

N 

^before = m% \2 V i + ( 20 '' 
i=l 

where (j>i — —GM/r t — ipi is the gravitational potential, and v CSCt i — \f2%pi is the escape speed for star i. The velocity kick 
relative to the center of mass is 

N 

Svi = Avi - Av = Avi - — ^ Av * ( 21 ) 

i=l 

where At; is computed for a continuous cluster in Appendix A. Star i is ejected if (vi + Svi) 2 > v 2 Bc i . The total amount of 
mass ejected from the cluster can then be formally written down as 

JV 

AM = -^mi0<(ej) (22) 

£=1 

where #i(ej) = 9 Uvi + Svi) 2 — f 2 sc ,i] is a useful shorthand notation; similarly 
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9i(b) = 8 [f^sci — i v i + ^^i) 2 ] f° r the bound stars. We neglect the possibility that stars unbound by our criterion remain near 
the cluster for a long time, possibly to become bound once again in a subsequent encounter with a black hole [Spitzer (1987) 
discusses orbits of this type] . 

To find the change in cluster energy we need the kinetic and potential energy changes, AT and AV, respectively. The 
change in kinetic energy arises both from the kicks to the bound stars and the loss of the energy of the ejected stars: 



JV JV 

AT = J-Y, " X ~ 

i=l i=l 



miVi6i(ej) + ^ rm (vi ■ Svi + ^vf^j 0»(b). (23) 



The change in the potential energy is AV = V a ft C r — V>cforc, where the potential energy before the collision is 

N JV 

V bcforc = -f V V ™* [fli(b)fl 3 -(b) + 2^(ej)^(b) + ^(ej)^(ej)] . (24) 

1=1 J^i=l 

and potential energy of the bound stars remaining afterwards is 



JV JV 

V aftcr = ~£ E ^(b)^(b)^ (25) 

i=l JV*=1 



for separation vectors fy. Consequently, 

JV JV 



A V 



2 

i=l i^i=l i=l j^i=l 

JV JV JV 



= -^9,(ejW.-jE £ ^N)^(ej)^» (27) 

j=l i=l 

where </3i still denotes the potential of star i from before the collision. The first term in eq. ( p7| ) is just the pre-collision potential 
energy of the ejected stars, and is 0(AM/M), the fractional mass loss. The second term is 0(AM/M) 2 , and will be smaller 
than the first as long as AM/M is small. In our calculations we neglect the smaller (and difficult to compute) second term, and 
evaluate the first using the King model potential at the positions of the ejected stars. To protect against gross inaccuracies, 
we terminate our simulations if | A M/M\ > 0.2 in a single collision. 
With these approximations, the energy change is 

N N 



AE = - V] TO: 

1=1 



(~«! + &(ej) + ^m, (vi ■ 5vi + ^8v^ ft(b). (28) 



It should be noted that the energy change due to ejected stars is always positive, and hence mass loss always heats the cluster. 
In the sum over bound stars, the v% ■ Svi term can have either sign and hence can heat or cool. The mean of this term is 
nonzero but small since it depends on the mass loss to make the final distribution slightly anisotropic; it also contributes 
substantial variance. The \8v* term gives rise to the familiar mean heating (King 1966; Binney & Tremaine 1987). 



5.1 The Model for Energy Input and Mass Loss 

For the analytic estimates and Fokker-Planck calculations of cluster survival, we shall need formulae for the mean energy and 
mass changes, AE and AM, their variances, and the cross-term AE A M. To get them we use a model for energy input and 
mass loss which is due to Chernoff et al. ( 198tj ) (hereafter called C.K.S.). In this model, the portion of phase space from which 



mass escapes is identified. The energy and mass of the remaining cluster are then found as integrals over the distribution 
function of the remaining stars. This model relies only on the information given by the distribution function from before the 
collision; no detailed solution of the Vlasov equation is attempted. The higher order effects which arise from the detailed 
alteration of the distribution function by the perturbing black holes are estimated and define the error in our method. 

We shall also extend the the C.K.S. model slightly by considering the fluctuations in the number of stars lost in any 
collision resulting from a finite number N of stars. Before the collision, the phase space density of the cluster is given by the 
King model in eq. (|l^). After the collision, we do not have a full expression for the phase space density including the effect of 
the change in velocity to all the stars, but we do know that all parts of phase space for which 8(e)) — 1 will have their mass 
ejected. Define /< = /km g #(b) and /> = /kin g #(ej) which are nonzero only for bound and unbound stars, respectively. The 
average of a quantity x over the bound, or ejected, stars will be written 

(x < = J r /' , (x > = \ / (29) 

/ d 6 r/< w / d 6 r/> v ' 
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respectively. If a symbol < or > is not specified, it means that the quantity is averaged over the entire phase space. 

The probability that a star is in the portion of phase space from which mass is ejected is p e j = (#(ej)}, which is an integral 
over the cluster model and is independent of N. We then suppose that the distribution for losing N e j stars out of the total N 
is given by 

P(iVcj) = ^0^exp(-p ej iV). (30) 
In averaging the expressions of we will only need the moments (N e j) = Np c j and (N^) = (iVpcj) 2 + iVpcj. 



5.2 Moments of AE and AM for Individual Collisions 



In this section,we integrate the first and second moments of AE and AM from §|| over the King cluster model and the 
distribution for the number of ejected stars P(A r c j). These averages will be denoted by an overbar to distinguish them from 
the averages over bound and ejected stars defined in the previous section. 
The mean mass loss and energy input are 



AM(&,??c,Vo) 



M 



-(%j)> 



(31) 



and 



AE{b,T) c ,i> ) 



\E\ 



M 

\E\ 



-Pcj 



~ W\ 

1 2 ^ 

2° + 



-P<y (-v 2 + <j>^ + (1 -p cj ) (v ■ Sv + t;8v 2 ^ 



+ ( v ■ Sv - 



l -Sv 2 



(32) 



Consistent with our neglect of terms oc (AM/M) 2 (and smaller), we only retain contributions at 0(rj 2 ) and (D(r]c)', note that 
p c j oc ijj, \6v\ oc T] c , and consequently (v ■ 5u)< - which would vanish in absence of mass loss - is oc 77^. To the same accuracy, 
the variance of AM is 



<rAM(b,r) c ,'4>o,N) _ p^ 
M 2 ~ N 

and the variance of AE is 

ai E {b,ric,4>o,N) M 2 



(33) 



E 2 



NE 2 



Sv + \sv 2 



- !>■ ; { [ \^ + < 



-Pej 



+ (1 — Pej) \V ■ Sv + -Sv 



M 2 



NE 2 

Finally, the mixed moment is 



((v-sv) 2 ) <+Pcj ^v 2 +<py^ 



(34) 



4AM(Mcio,iV) _ (AM A E)(r] c , 4>o, N) AM{r} c ,i> ) AE(rj c ,iPo) 



M\E\ 



M\E\ 



M 



\E\ 



M 

7j\e\ 



(35) 



keeping only terms 0(rfc) and larger. 

We have calculated the integrals needed to evaluate these moments u sing Monte Carlo methods. The velocity kick is 
evaluated using eq. (^) and the center of mass velocity kick is given in eq. (/II). A sample of the results for ipo = 2.0, 4.0, 6.0 
and rj c = 10~ 2 are plotted in Figs. |^, |^, and ^. For ?/>o = 2.0, 4.0, and 6.0, the core radii are at r cola /r t = 0.2, 0.1, and 
0.05, respectively, while the half mass radii are at rum/rt = 0.3, 0.2, and 0.15, respectively. The Monte-Carlo integrals were 
computed to an estimated fractional error of 5 per cent. 

The separate contributions to the energy input in eq. jj^ have been plotted in Fig. [| The contribution of mass loss is 
comparable to that due to the Sv 2 /2 heating usually considered. Note too that the (v ■ Sv) term is surprisingly large for small 
b encounters, contributing ~ 5 per cent to the energy input for impacts near the cluster center. The variance in eq. (^) is 



also broken up into the heating part, ((v ■ Sv) 2 ^ and the mass loss part, p e j ( (f 1,2 + 



which are comparable in size. 
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5.3 The Tidal Limit 



The calculation of AE, AM — vNctam, "asi an d °"aeam i s simplified considerably in the tidal limit, partly because 



mass loss is restricted to particles with speeds very close to the escape speed. Spitzer (1958) computed AE without mass 
loss, however, mass loss contributes significatly to AE, sometimes exceeding the Spitzer term. The expressions found below 
illustrate the importance of mass loss explicitly and are also useful for understanding disruption of clusters by high mass black 
holes. 

Formally, we should expand in the two parameters rj c and rt/b <C 1. Here, we keep only the lowest powers of r] c and 
j3 — rt/b. In this approximation, we find 



M = ~ Cm W (36) 

where 

Cm = /^d^e^) 372 , (37) 



9^C t 4 tec 
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io- z o.i 1 

b/r t 

Figure 3. |(AM(d,rjc))|/M. 



and £ = r/a, & = r t /o, a = (* /47rGp ) 1/2 , jtit = M/M , M = iirpoa 3 , f c = p /pi, g c = ffg/fff, and 6(0 = *(£)/*(> [pi 
and cti are defined in eq. @]. Hence a 2 AM /M 2 = p e j/-/V = (Cmm /N)(t] 2 /f3 4 ') where CWm = CW- To 0(77?), 



AE(b, Vc ,$ ) Mr GM II 2 \ X GM 2 4 G 2 M 2 h M{r 2 ) _ ^ 

MrT (2 // Pcj ^ + 3 TO&i = c ^ (38) 



i^i 

with 



GM 2 4 G 2 M 3 (r 2 ) 

Ce = Cm — r^j- + — — r=q — 3 — • (39) 



r t |£| 3 \E\a 2 r 4 



In addition to the Spitzer energy loss {5v 2 /2), Cb includes the energy carried away by the massloss, AM = — p e jM, from near 
the escape surface, where the energy per mass is —GM/rt- Note that the ejected mass can come from any spatial position 
in the cluster as long as the velocity is close enough to the escape velocity. Indeed, the mean radius from which mass is lost 
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Figure 4. N<r^ E (b, rj c )/E 2 . 



depends weakly on cluster concentration, and is roughly 2r t /3 for < ipo < 8.5. The variance to the energy input is 



<?A, B (Mc,V'0, AO 



E 2 

with 

Gee = Cm 



( GM 2 \ , 8 1 G 2 M 2 h M 2 . 2 a 

(r v ) 



A^' J V |S|r t 



+ 



9 AT S 2 V; 2 el 6 4 



N f3 4 



(GM 2 



8 G 2 M 4 , 2 2X 



V|S|r t y ' 9£ 2 a 2 r t 4 
and the mixed moment is 

^asam(Mc, Vto, AT) M 



M\E\ 



N\E\ 



M 



GM 



~N\E\ Pci n 



Gem ijc 
N ]F 



(40) 



(41) 



(42) 



with Gem = C M {GM 2 )/(\E\n). 

Another extremely important quantity will be the change in v = (Er t )/(GM 2 ) in the tidal limit. Using the fact that the 
cluster is tidally limited, v can be expressed as v = (r t , /(GM 1/3 )) x (E/M 5/s ) where the subscript refers to some reference 
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Figure 5. Dimcnsionless coefficients for fractional energy input and mass loss in the tidal limit. 



value. The change in v is then 

' AE 5 AM" 
-- 



A v 



M ( 



M [c e --Cm )n 



E 3 M J r 1 V"" 3' 
which can be written in the usual form 



Av(b,r) c ,ip ) 



_ n '/c 



with C v = Ce — 5Cm/3. The variance of the mean change in v is given by 

cL,(Mc,i/>o,iV) 



2 

°~AE 
E 2 



2 2 

5 caeaa/ 25 Q- am _ Ci_ 
3 9 A/ 2 ~ N 



t'V f] c 



(43) 



(44) 



(45) 



where C vv = C E e - 5C E m/3 + 25C M m/9. 

The seven coefficients for the tidal limit are plotted in Figs. |H| and|^ as a function of ipo. All the curves show a monotonic 
decrease as ipo increases over the range Vo G (0.0, 8.5). Most noticeable is that C v becomes less than zero at roughly ipo — 5.5. 
Since v = v(ipo), whether an impact will make the cluster more or less concentrated is determined by the value of ipo for that 
cluster. Clusters with ipo < 5.5 are driven toward dissolution and clusters with ipo > 5.5 are driven toward core collapse. This 
agrees well with CKS who derived the same result in the context of perturbations by giant molecular clouds. 

For penetrating encounters, the situation is quite different. Over almost the entire range of ipo and r/co, the diffusion 
coefficients yield positive changes in v, implying dissolution. 



5.4 Evaluation of the Diffusion Coefficients 



In §5.2, we found AE, AM, o\ E , o~\ M , and <j\ E am f° r individual collisions as a function of impact parameter b and relative 
velocity V le \. When Aibh is small enough that the changes in the cluster properties are always small for single impacts, then 
we may average over many encounters to find the changes in E and M over a time period which includes many collisions, but 
for which the changes in cluster properties are still small. 
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0.5 - 



Figure 6. Dimensionless coefficients for fractional changes in u in the tidal limit. 



We weight the formulae for the various moments by the differential rate of collisions 

dr = n bh 27rta&/(V r el)14eldV r el 

and integrate over b and V xe i where /(Kci) is given in eq. ([ji]). The results may be written 
D E = dT'AE = T \E\d E 



(46) 



D M = / dTAM = -T Md M 



D EE ■■ 
Dmm = 
Dem = 



dFa&E = jrjT E 2 d E E 
dr<ji M = — T M 2 d M 
dra| BAM = -^T \E\Md EM 



(47) 



where To = nbh^r^Vc, and d_e, dM, dEE, and dsM are dimensionless functions of ipo and r] c o = (Mbh/M)(ao/V c ). Note that 
minus signs have been inserted in Dm and Dem to make dM and dsM positive. 

The diffusion coefficients d_B, d*/, dEE, and dsM involve eight-dimensional integrals over phase space, b, and V re \. Inte- 
gration over two velocity angles is elementary, and the integral over velocity magnitude was done using Runge-Kutta. The 
remaining five integrals (three spatial coordinates, b, and V le i) were done via Monte Carlo (to an accuracy of 5 per cent). The 
value of 6 max was taken to be 10r t . 

For use in our Fokker-Planck calculations, the four diffusion coefficients were tabulated for almost all the values in the 
range tpo — 0.0 — > 8.5 and 7/ c o = 0.0001 — > 0.5. Due to the extremely long integration times for small values of rjco and 
large ipo, a few values were obtained using analytic fitting formulae of the form = ke ln(l/?7co)??c(b d-M = km ln(l/r] c o)r/ 2 , 
dEE = k,ee ln(l/?7 c o)?7co, an d &em = hem \n{l/rj c o)ri 2 , which represented the data well; the coefficients for these fits are 
given in Table 2. The fits were generally good to <5 — 30 percent and are useful for quick estimates, although for more precise 
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Table 2. Fitting Formula Coefficients for d^, d^, d^Ej an d ^EM ■ 



$0 


k e 


k E 


km 


km 


kee 


K-EE 


KEM 


K>EM 



0.0 


12.8 


70.5 


4.14 


22.8 


20.8 


114 


7.29 


40.1 


0.5 


13.4 


80.5 


4.34 


26.1 


23.3 


140 


7.90 


47.5 


1.0 


11.0 


73.2 


3.38 


22.5 


18.3 


122 


5.89 


39.2 


1.5 


10.1 


75.3 


3.10 


23.1 


16.8 


126 


5.31 


39.6 


2.0 


7.05 


60.1 


2.06 


17.6 


10.0 


85.3 


3.07 


26.2 


2.5 


7.68 


76.2 


2.24 


22.3 


12.2 


122 


3.56 


35.4 


3.0 


6.42 


75.8 


1.85 


21.9 


9.64 


114 


2.72 


32.2 


3.5 


5.42 


78.2 


1.51 


21.8 


7.95 


115 


2.09 


30.1 


4.0 


4.85 


88.0 


1.39 


25.2 


7.33 


133 


1.93 


35.0 


4.5 


3.80 


89.5 


1.06 


24.9 


5.55 


131 


1.33 


31.3 


5.0 


2.85 


90.2 


0.740 


23.4 


4.05 


128 


0.985 


31.2 


5.5 


2.33 


102 


0.702 


30.8 


3.36 


147 


0.780 


34.1 


6.0 


1.88 


117 


0.603 


37.7 


2.73 


171 


0.631 


39.4 


6.5 


1.40 


126 


0.499 


45.0 


2.03 


183 


0.489 


44.0 


7.0 


1.06 


137 


0.379 


48.8 


1.56 


201 


0.329 


42.4 


7.5 


0.836 


141 


0.363 


61.4 


1.24 


209 


0.308 


52.0 


8.0 


0.677 


131 


0.320 


63.2 


0.971 


192 


0.243 


48.0 


8.5 


0.647 


126 


0.362 


70.7 


0.928 


181 


0.277 


54.0 



Note: The dimensionless diffusion coefficients are approximated by d = ftln(l/?7 c o)J7j?Q and d = 
k\n(l/i) c o)fil where d = d(r t 2 /rg m ) and fi c0 = r) c0 (GM /r^a^) 1 / 2 . 

numerical work interpolation on tabulated values was used. The values generated using the fitting formalae were used only for 
the very smallest values of r/co for the cluster PAL 5, where the fitting formulae are most accurate. Fig. ^displays the various 
coefficients for several values of ffco which are shown increasing from bottom to top. The triangles represent data which was 
generated by the Monte Carlo program itself, while the open circles represent the data estimated using the fitting formula. 

It can be seen that for a fixed value of rj c o, the diffusion coefficients in Fig. |^ vary by about an order of magnitude over 
the range — 0.0 — > 8.5. The coefficients for the fitting formulae show a similar range of variation. This spread can be 
reduced considerably if we recall that many quantities vary little over the King sequence when expressed in terms of the 
half-mass radius, rh m - Choose a new unit of rate fo = n^nr^Vc and a new collision parameter f/co = ilco(GM /r^Po) 1 2 . 
The constants k in the new dimensionless diffusion coefficients d = khi{l / i) c o)i)% = d(rt/r 2 m ) then vary by about a factor of 
two over the range ipo = 0.0 — * 8.5, as can be seen in Table 2. 



6 'SLOW HEATING' LIFETIMES 
6.1 The Fixed Cluster Approximation 

When Afbh << Mhigh, the properties of a given cluster only change slightly over many encounters with black holes. In this 
limit, cluster survival implies an upper bound to /haio-Mbh (up to logarithmic corrections), or Mhh assuming a certain value of 
/halo- Very rough estimates of the critical black hole mass for survival over time T can be made if we neglect the evolution of 
the cluster profile over time T, and require the mean energy or mass change not to exceed some threshold over that time span 
(e.g. Moore 1993). The values found here are intended as benchmarks for comparison to the more detailed results, obtained 
in the next section, that include cluster evolution and statistical spreads in the energy and mass changes in time T. 

For a fixed cluster profile, dE/dt, dM/dt, and du/dt are independent of time, and {AE(T)} = Fo\E\d E T and (AM(T)) = 



-ToMdhiT , with the coefficients and dM defined and computed in §5.4. Choosing a time of T = 10 years as the time over 



which the clusters have been subjected to collisions, we defined three distinct values of Mbh.crit by {AE(T, Mbh.crit)) = \E\/2, 
{AM (T, Mbh.crit)) = -M/2, and {Av{T, M bh ,crit)} = [D E /\E\ - (5/3)(D M /M)]\v\T = -0.6 - v. The results are presented in 
Table 3 for the sample of nine weakly bound clusters employed by Moore (1993). Moore's results correspond most closely to 
Mbb,ait(E) as he used the criterion {AE(T = 7 x 10 9 years, A/bh.crit)) = 1^1- When comparing to Moore's results one should 
multiply Mbh,crit by a factor of 2/0.7 to account for the larger value of T and smaller AE used here. Except for PAL 13, our 
results agree with Moore's to within a factor of order a few, which is resonably close considering the improvements made here 
(e.g. inclusion of mass loss, correct V Te i distribution). 
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Figure 7. Dimensionless diffusion coefficients. The triangles represent Monte Carlo results. The circles represent the incomplete data 
given by the fitting formulae. The values of T) c o shown are 0.0001, 0.00025, 0.0005, 0.00075, 0.001, 0.0025, 0.005, 0.0075, 0.01, 0.025, 0.05, 
0.075, 0.1, 0.25, 0.5, increasing from bottom to top. 



Table 3. Cluster Destruction for a Fixed Profile 



Cluster Name 


M h ^ clit (E)/M Q 


M bh , crit (M)/M 


M bhiCrit H/M 


AM 4 


1500 


7000 


710 


ARP 2 


3100 


13,000 


4900 


NGC 5053 


4100 


18,000 


5200 


NGC 7492 


3400 


15,000 


6500 


PAL 4 


170,000 


8,500,000 


350,000 


PAL 5 


1100 


4300 


1100 


PAL 13 


3100 


13,000 


6300 


PAL 14 


14,000 


68,000 


15,000 


PAL 15 


9100 


39,000 


6100 



Note: Afbh crit(B)/Af0 is the critical black hole mass for the criterion SE/\E\ = 0.5. 
A^bh cT\t(E)/MQ is the critical black hole mass for the criterion \SM\/M = 0.5. M^ cy \ % (v)/Mq 
is critical black hole mass for v to go from its initial value to v = —0.6. 



6.2 Gaussian Model with Evolution 

In jj js.lj two approximations were made. First, the cluster was treated as having a fixed profile for which the diffusion coefficients 
did not change over time. The second approximation was that the energy input and mass loss had sharply defined values over 
any time interval. In this section, the evolution of the cluster is followed over appropriately chosen intervals St and a Gaussian 
distribution of energy input and mass loss is assumed. 
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The characteristic time for iV m i n collisions to occur inside 6 max is 

T coli = ^2™, „ oc M bh R 2 g N min . (48) 

As most of the energy input is given by the penetrating encounters, a conservative estimate is to set 6 m ax = rt in T co n. If e max 
is the largest cha nge allowed for the cluster in a time period, then a time Tchange over which the cluster evolves significantly 
is given by (§5.4) 

R 2 

Tchange = 2 OC . (49) 

The Gaussian approach is justified if there exists a St such that T co ii < St < STchmge, so the cluster properties change little over 
St but enough collisions occur that the distribution of energy inputs and mass loss is approximately Gaussian. An estimate 
of the largest black hole mass, Mbh,f P , for which the Gaussian approach is justified is 

M bMp < M^(^X' 2 , (50) 

or VcO ?S (e m ax/A r m i n ) 1//2 . The critical mass can only be found in the present 'diffusion' approximation when e ma x <C 1 and 
iVmin > 1 so that the cluster must be destroyed for 7/ c o -C 1. 

Let the vector r = (n,T2) = (AE, AM) where AE and AM are the energy input and mass loss over the time interval St; 
the expected values of AE and AM over the time St are (ti) = (AE(St)) and (t 2 ) = {AM (St)} . The normalized distribution 
of energy input and mass loss is then 



P(T) d 2 r = V^gg d a Tatp |_l (r _ {T) f S{T _ (r)) | , (51) 

where 3(5t) is the covariance matrix. The matrix H _1 may be shown to be 
((n-(n)) 2 ) ((n - (7i»(7i - <r 2 ») 



((ri - (rr))^ - (r 2 ») ((r 2 - (r 2 » 2 ) 

inverting implies 

1 f <(t2-<t 2 » 2 > -((n - (n))(r 2 - <r 2 ))) 



" det^" 1 ) V -((n - (n»(r 2 - (r 2 »> ((n - (n)) 2 ) 

These expressions reduce to the standard sum of Gaussian terms when the correlation ((n — (ti))(t 2 — fa))) is zero. The 



eigenvalues of S are A± = \ Tr(3) ± yj (Tr(H)) 2 -4det(S) . 

In order to choose values for AE and AM, we define a new set of variables £ = (£+,£_) along the eigenvectors of S; the 
distribution of £ is 



P( C )d 2 C = ^±Ed 2 Cexp{-A + a/2-A_CV2}, (52) 
and t — (t) = A£, where the unitary matrix 

A i+ Ai_ \ 1 / Hi2 A- - H 22 

A 2+ A3- J ~ y /(\_ - =22) 2 + H 2 2 V "( A - - H 22) 

After choosing (± from this Gaussian distribution, we multiply by A to get 
AE = (AE) + Ai+C+ + Ai_C- 

AM = (AM) + A 2+ C+ + A 2 _C-- (53) 

In the limit of infinite numbers of collisions N c in time St (the limit in which the Gaussian model is equivalent to a 2-D 
Fokker- Planck equation) 

(n) = r \E\d E St 
(r 2 ) = -r Md M St 

Su(5*) = ^T E 2 d BB 8t 
^2 (St) = ~r \E\Md EM 5t 

=■22 (St) = lr M 2 d M M5t. (54) 
For finite N c , these expressions are accurate to 0(N^T 1 ). 
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Let us ignore the correlations of the change in mass and energy and focus solely on changes in cluster energy in order 
to get a qualitative feel for the 'distribution of cluster states'. The Fokker-Planck equation for the distribution of energies, 
P(E,t), of an ensemble of clusters is 



If wc ignore changes in the cluster properties and take De and Dee to be constant, this equation has the solution 



P(E,t) = 



(2ttZW)i/2 



exp 



(E- 



D E tf 



2D EE t 



(55) 



(56) 



where we have chosen the initial condition E — Eq at t — 0. For small times t < D EE /D E , the width of the Gaussian is 
much larger than the mean implying a very spread out set of cluster states, and visa versa for times t > D EE /De- As we 
have already shown in §2.3, for times t ~ Tdismpt the set of final states is sharply defined about the mean. Using the survival 
criterion AE < f\Eo\, wc find that the probability of survival after time t is given by 



Ps = 



E +f\E \ 



dEP(E,t) = - 



1 + erf 



f\E\-D E t\ 



(57) 



V2D EE t 

The same equation would have resulted if we had used the two-dimensional Fokker-Planck equation for E and M. In terms 



of the fitting formulas of §5.4, the characteristic mass range, <5Mbh,, 

SMbh,cvit 



over which P s changes from one to zero is 



K-EE 



< 1 



(58) 



M b h,crit V NfK E 

so that in the Mbh <C Mhigh limit the value of Mbh,crit found from just one history is likely to be close to the statistical mean. 



6.3 Simulations and Results 

We chose the accuracy parameters from the previous section to be e max = 0.01 and iV m m = 10. The time step was chosen 
so that the estimated change in energy would be 0.01. As the diffusion coefficients tended to increase as ipo — > and M 
decreased, the time step was decreased if two events in a row occured with fractional changes in v, M, or E greater than 0.01. 
Also, the time step was increased if the estimated number of events in St was less than N m i B . If Mhh > Mbh,f p , as described 
in the previous section, the simulation was stopped because the approximations had broken down. 



For each of the nine weakly bound clusters from Moore (1993), values of Mbh were chosen and the cluster was evolved 
to one of the following outcomes: (1) the Fokker-Planck assumption broke down; (2) T = 10 10 years was reached and the 
cluster survived; (3) the quantity v went out of the range (—2.13, —0.6), signalling either core collapse or 'dissolution'. Survival 
probabilities and their uncertainties were derived from the number of surviving clusters out of 1000 simulations using Bayesian 
arguments. The results are shown in Fig. ^| for Moore's clusters as well as one more 'normal' cluster. Note that all curves 
asymptote to /dost = 1-0 as Mbh — > oo because single-event destruction has been ignored in these simulations. 

PAL 4 did not satisfy the Fokker-Planck assumption at the black hole mass required to disrupt it in T = 10 10 years. Since 
PAL 4 lies at such a great galactocentric distance, encounters are so infrequent that disruption in 10 10 years requires individ- 
ually destructive collisions. Hence a critical black hole mass Mbh.crit cannot be determined by the slow heating approximation 
for PAL 4; the more detailed collision-by-collision encounter history of §^ is warranted for this case. 

As for the rest of Moore's clusters, we see that their values of Mbh,crit agree well with the last column of Table 3; the 
values obtained from the Gaussian Monte-Carlo simulation are lower by less than a factor of two. As clusters are heated and 
lose mass, diffusion coefficients increase as ipQ — > 0, and rj c o decreases leading to smaller Mbh, cr it than for fixed cluster structure 
and mass. Also note that the fractional width, <5Mbh,crit/Mbh,crit, of the 'transition region' over which /d cs t = 0.1 — » 0.9 varies 



as TV 1//2 , as predicted in §2.3 



We also simulated a large cluster with M = 1O 6 M0, TV = 1.4 x 10 6 , r t = 50 pc, r CO rc = 5 pc, and R g — 10 kpc; hence 
rjjo = 4.85 and the central one-dimensional velocity dispersion ao = 9.64 km s _1 . Fig. ^ shows P s as a function of Mbh for this 
cluster. The critical black hole mass is quite a bit larger than that found with the small clusters from Moore (1993) because 
M is larger. Note that Mbh.crit ~ M/10 for this cluster, which is on the verge of the limit in which many black holes will be 
inside the cluster at any given time. 



Klesson & Burkert (1995) recently argued that Moore's determination of Mbh.crit for his set of nine diffuse clusters was 
flawed because deviations from the mean energy input for each encounter history would be quite large, implying a significant 
scatter in the values of Mbh.crit- They contend that clusters are destroyed as a consequence of a small number of encounters with 
large energy inputs. Instead, we find that for eight of Moore's nine clusters, all except PAL 4, disruption is due to numerous 
encounters with b < r CO re which cause small changes individually. For these light clusters, the expected number of impacts, 
TV core (Afbh.crit), inside r core is large for Mbh = Mbh, cr it- We find that five of Moore's clusters have A r core (Mbh,crit) > 200, 
two have A r core (Mbh. C rit) > 50, one has /V CO rc(Mbh.crit) = 20, and PAL 4 has A r corG (Mbh,crit) ^ 1. So (aside from PAL 4) the 
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Figure 8. Fraction of clusters destroyed in T = 10 10 years plotted against black hole mass. These curves were generated using the 
Gaussian Monte Carlo method. 



expected fluctuations in the A r corG (Mbh,crit) are <20 per cent for eight of Moore's clusters, and <10 per cent for seven of them. 
Only for PAL 4, which sits at a large galactocentric radius where few encounters occur, is the evolution near Afbh,crit very 
stochastic. This case will be discussed in the next section. 



7 THE LARGE M bh LIMIT FOR THE SURVIVAL PROBABILITY 
7.1 Theory 

The previous section treated the small A/bh regime in which many collisions slowly add energy and remove mass from a cluster. 
This approximation will be valid for M b h « Afhigh, the mass at which the cluster can be disrupted in a single collision. In 
the opposite limit for which Mbh ^ Afhigh, any collision within a 'destructive impact parameter', bd, will destroy the cluster. 
In addition, the cumulative effect of many collisions outside 6d can also destroy the cluster. 

Formulae for the changes in E, A/, and v for a single collision in the tidal limit were derived in §5^; see eq. (0), 



and @. The three formulae have the same scaling with r\ a and b and only differ by a factor which depends on ipo. Hence, 
in the following derivations we only use the formula for energy input to derive results explicitly. To find the expressions for 
disruption by mass loss or change in v one must only substitute C'e — * Cm, C v . 
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The energy input can be rewritten as 

AE - 4 ' *- - ^ 4 



\E\ 



where the impact parameter for AE/\E\ = / is 

1/4 



fed 



1/2 



Ce 

f 



(59) 



(60) 



This energy input is sharply peaked about the mean for fe ~ fed because the impact parameter at which AE ~ <tae is larger 
than fed by a factor of TV 1 / 2 . 

A cluster has been destroyed if an impact occurs with 6 < fed. The critical impact parameter fed > r t for Mbh > Mhigh = 
(MV rcl /a )(f/C E ) 1/2 - Using (V rel ) ~ 1.47V C [see eq. @], and (r ~l - 10km s" 1 implies M high ~ 10 - 100M. 

The mean number of destructive encounters, Nd, which occur in a time T is given by the integrated rate of encounters 
inside b d (V lel ): 



N d (T) = im hh T J d^. ol /(Kci)Kcib d (Kci) = (^y) ^r 2 t a T. 



(61) 



Note that Nd{T) oc /halo but is independent of Mbh- Consequently, when Mbh > A/high, limits can be derived for /halo but 



not M bh (Wielen 



For given / and T, the Poisson probability of no individually destructive events is given by 

Ps (/halo) = exp[-JV d (/ ha l )]. 



(62) 



The expression for P s becomes more complicated when we include the cumulative effect of all the nondestructive encounters 
outside fed. As most collisions have 6 >> fed, they do not significantly change the cluster properties. These gentle collisions 
can be modelled using diffusion coefficients. 

The mean rate of change of E due to impacts outside fed is given by 

fed(Kd)" 



E 

W\ 



dV ro i/(Kci)Kei 



f- 



6 4 



27rfedfe 

where Td = Nd(t)/t. Similarly, the time derivative of the variance is given by 



E 2 



dV re l/(Kel)Kel 



/ 27rfedfe 


" f 2 &d (Vrel)" 
1 b 8 







= \fv d 



(63) 



(64) 



where we have ignored the small variance per collision; this variance is due only to the variable number of encounters in 
a given time period. In addition, collisions inside fed act as a 'sink' for clusters. The appropriate equation for P(E,t), the 
probability density at time t and energy E, is 



dP _ p9P 1 . 2 d 2 P 

— - -L d P-E— + -a E — 

where again the diffusion coefficients are treated as constants for simplicity. The solution with E = Eq at t = is 



P(E,t) 



exp[-T d t] 
(27ra|t)V2 



exp 



(E - E - Etf 



2&%t 



(65) 



(66) 



This expression reduces to the Mbh < Mhi g h solution when fed 
probability that the change in E has been less than f\E\, or 



1 

2 6 



1 + erf 



f\E\ - ET 



= 2 6 



-N d 




0. The probability of survival is then the cumulative 



(67) 



-1/2 



exp(— 5Nd/2), which decreases much faster. 



For small Nd, P s — exp(-Nd), as in eq. (p2|), but for large Nd, Ps - 
Analagous results can be derived for changes in v and M. 

This decrease in the survival probability bel ow th e simple exponential decline due to the effects of distant encounters has 
been noted by previous investigators. In Wielen ( 1988 ), Fig. 5 presented a numerical determination of the survival probability 
in the large Mbh limit for clusters bombarded by both black holes and giant molecular clouds. Using Nd(T) ~ 0.5 x T/7i/ 2 , 
where T x / 2 is the time over which half the clusters have dissolved, ou r eq. ( |67j) is a rather good fit to the figure. 

The results of this section do not agree with Klesson & Burkert ( 1995 ) since the curves of P 3 against Mbh in their Figs. 
8, 9, and 10 show no systematic evidenc e of a n asymptotic P 3 independent of Mbh. This is surprising as they have included 
exponential adiabatic damping [Spitzer ( 1958 ), but see also Murali & Arras (1998) and Weinberg (1994) for a more recent 
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Table 4. Limits on /halo m the Large M D h Limit 



Cluster Name 


Criterion^ ) 


M high /M 


/haloes = 0.9) 


/halo {Pa = 0.5) 


/haloes = 0.1) 


AM 4 


V 


3.8 x 10 5 


0.028 


0.14 


0.33 


ARP 2 


V 


1.6 x 10 7 


0.014 


0.07 


0.16 


NGC 5053 


V 


1.8 x 10 7 


0.011 


0.06 


0.13 


NGC 7492 


1/ 


1.5 x 10 7 


0.029 


0.15 


0.34 


PAL 4 


1/ 


1.4 x 10 7 


0.174 


0.90 


>1.0 


PAL 5 


V 


8.6 x 10 6 


0.005 


0.02 


0.05 


PAL 13 


V 


6.4 x 10 B 


0.053 


0.27 


0.63 


PAL 14 


V 


8.7 x 10 6 


0.035 


0.18 


0.41 


PAL 15 


V 


4.4 x 10 6 


0.025 


0.13 


0.30 


AM 4 


20%M 


4.1 x 10 5 


0.031 


0.20 


0.67 


ARP 2 


20%M 


7.6 x 10 6 


0.006 


0.04 


0.14 


NGC 5053 


20%M 


1.0 x 10 7 


0.006 


0.04 


0.13 


NGC 7492 


20%M 


4.9 x 10 6 


0.009 


0.06 


0.20 


PAL 4 


20%M 


8.9 x 10 6 


0.109 


0.72 


>1.0 


PAL 5 


20%M 


5.8 x 10 B 


0.003 


0.02 


0.07 


PAL 13 


20%M 


2.1 x 10 6 


0.018 


0.12 


0.39 


PAL 14 


20%M 


6.1 x 10 6 


0.025 


0.16 


0.54 


PAL 15 


20%M 


3.9 x 10 6 


0.023 


0.15 


0.50 



(a) The criterion 'j/' means that neither single destructive events nor many nondestructive 
events caused v —0.6. The criterion '20%M' means that no single events of 20 per cent mass 
loss occured. 



version of these effects] which would decrease the energy input per encounter and hence increase P a . The lack of encounters 
with small energy input in Klesson & Burkert's ( 1995] ) simulations could result if their method of choosing the collision 
parameters oversampled small b and V le i compared to our method. 

Table 4 contains the allowed values of /halo derived from the fixed cluster approximation of this section for the nine 
clusters found in Moore (1993). Two different destruction criteria were used: (1) v out of bounds in T — 10 10 years due to 
the combined influence of single destructive collisions and multiple nondestrucive collisions [eq. (^57j)], and (2) a single episode 
of 20 per cent mass loss in T = 10 10 years [eq. flffi^)]. Also given are the values of Mhi g h, the black hole mass above which 
the cluster can be destroyed in a single collision, appropriate to each criterion. Three values were used for the probability of 
survival: P s = 0.1,0.5, and 0.9. 

For 50 per cent survival probability using the criterion v — > —0.6, the limiting values of /halo range from 0.02 (PAL 5) to 
0.90 (PAL 4). It is unlikely that /halo > 0.3 since P s < 0.1 for most of the clusters in that case. These limits depend sensitively 
on galactocentric radius and cluster size, as is evident in the scatter in the values of /halo- 

The results of this section will be tested in §[] using the results of the full Monte Carlo simulations. One wrinkle which 
appears in the results is that the different criteria for destruction can compete with each other so that the results of this section 
are not always good approximations to P 3 (Nd). Only in the cases where one method of destruction completely dominates over 
all the others do the results agree accurately. 



8 MONTE CARLO SIMULATIONS OF INDIVIDUAL COLLISIONS 
8.1 The Set-up 

Previously, Mbh.crit and /halo were calculated for the two limiting cases of slow heating (small Mbh) and single-event destruction 
(large Mbh). The major simplification for both methods was the neglect of detailed, collision-by-collision evolution of the 
cluster. Both calculations relied on approximations. The Gaussian Monte-Carlo method relied on restricting collisions to very 
small energy and mass changes individually, and involved averaging over many collisions to streamline the computations. In 
this section we present a Monte-Carlo simulation of individual black hole encounters. We make no approximations for the 
energy input besides the impulse and straight line approximations; consistent with the former, we neglect displacement of 
cluster stars during the encounter. This treatment is simpler than an N-body simulation because the evolution of the cluster 
is mapped by the sequence of King models. As before, the King sequence limits the evolution we allow. 

The calculations amount to simulating the 'Green's function' of the cluster. A cluster in a given initial configuration is 
subjected to collisions with the halo black hole population. After each collision the King model for the cluster is altered, using 
its post-collision E and M, and assuming r t oc M 1 ^ 3 . An evolutionary history of the cluster is mapped out over 10 10 years, 
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and a final cluster state is found. As this process depends on many random variables, a range of final states is possible for 
each initial state, and many realizations are needed to find the Green's function. 



8.2 The Probability Distributions for t, b, and V le \ 

In we derived the distribution of relative cluster-black hole speeds assuming an isothermal black hole halo. For the Monte 
Carlo simulations, we need to know the probability that the next collision suffered by a cluster occurs a time interval t to 
t + dt after a given encounter, and involves an impact parameter in (b, b + db) and a relative speed in (V Te i, V Te \ + dVrd). Since 
the number of collisions inside b is oc b 2 , we chose a maximum impact parameter b max , and define <r ma x = Tr&max- Let us first 
consider collisions with a single relative velocity, V Ie i- The expected number of collisions in a time t is 



N 



Kddr, 



so the Poisson probability that there are no collisions for the time interval (0, t) and then one collision in (t, t + dt) is 



dPi(t) = exp 



Kddr 



Vreldt. 



Multiplying by the probability 2nbdb/a nla _ x that the collisions is in (b, b + db) gives 

f 2ivbdb 1 



dPi(t,fc) = < exp 



7ibh0"maxKcldr 



Vteldt > X 



I 0" n 



(68) 



(69) 



(70) 



The first bracketed factor can be used to choose t, while the second bracketed factor can be used to choose b. 

Next, suppose there is a set of black hole populations, each with a single relative velocity V Ic i,i', let fi be the fraction of 
all black holes with V te i^. Then the probability that the next collision has relative speed V re i : i is 



dPi (*,«,&) 



exp 



^bh^max (V rel )dT 



,. J fiVr C l,i \ w ( 27Ybdb \ 

"-bhO- m ax(Kcl)dt } X <^ — } X < } 

\Viel) I I cr max J 



(71) 



where (V Te i) — fiV rc \ y i. Take the continuum limit by substituting fi — > f(V re i)dV re i to get 



dPl(t,6,Kel) 



exp 



(Kel)dT 



^bh^max 



2-Kbdb 1 



f /(Vrol)VroldKcl 

L (T max ) 



(72) 



The values of t, V Te i, and 6 may be chosen from the distributions in the first, second, and third bracketed factors, respectively. 

Note that our derivation does not require T = nbhO"max(Koi) to be time independent, so eq. ([72]) holds for any orbit, not 
just circular ones. Also, the velocity distribution found in the middle brackets is not f(V Te \), but is the distribution of speeds 
measured for the first collision after some chosen reference time. For example, the mean of this distribution is 



3/2 + 2: 



(73) 



( x + s) erf ( :E ) + 7¥ ex P ^ X ' 2 

where x — V c i/V c - As x — > 0, V IC \/V C —* 3v^"/4, and as x — » 00, V rc i/V c — > 1; for circular orbits eq. ( |l7| ) and eq. ( [7^ ) imply 
V^2{x = l)/(V le i(x = 1)) ~ 1.16 , or, V^i{x = 1) ~ 1.70 x V c . 

For globular clusters on circular orbits, the rate of encounters, T = nbhO"max(V rc i) inside a- max is constant. The normalized 
cumulative distribution for the collision time is then P c df,t(i) = 1 — e~ Tt . 



8.3 Evolving the Cluster 

Distinct histories are computed for individual clusters at fixed R g with identical initial values of M, r t , r COTC , and N, the 
number of cluster stars. The mass per star is held fixed at m = (M/N)\t=o even as M (and N) evolve; multiple species with 
different particle masses and spatial distributions were not considered. The maximum impact parameter & ma x is also held 
fixed. 

At any time t > 0, the cluster has known values of M, r t , and E which allow determination of ipo, and hence the smooth 
King model phase space distribution. From this distribution, we choose random positions and velocities for the stars in 
the cluster at the time of its next encounter with a black hole, t ncxt > t (to keep fluctuations to a minimum, we actually 
choose velocities first in symmetric pairs ±i> and then choose positions from their distribution given the velocity) . The values 
of tnoxt, 6, and V le i are chosen according to eq. (fr|), and velocity kicks Av t are computed for each cluster star using eq. [j| 
(by definition b = be x and V rc i = V vc \e z ). 



© 0000 RAS, MNRAS 000, 



24 P. Arras and I. Wasserman 

Given the Avi, the center of mass velocity kick is Av = (1/N) Avi\ even stars which are ejected are included in 

the sum. Star i is ejected from the cluster if (vi + Ad, — Aw) 2 > 2ip(ri), where V( r is the pre-collision potential at r». If 
JV' stars remain in the cluster then M' = M(JV'/JV) is the new mass of the cluster. The new energy of the cluster is the total 
energy input to the remaining stars (i = 1, 2, JV') plus the contribution from the ejected stars (j — 1, 2, TV — JV'); 



The new parameters r[ and ip' are then found by the steps outlined in §|3j. In this way, we have a new nondimensional model 
of the cluster after the collision. 

These steps are repeated until either t = 10 10 years or the cluster 'dies'. We keep track of six different conditions for 
destruction. For two of these, the simulation is stopped at once. The first, l v out of bounds', triggers if v goes out of the 
range (—2.13, —0.6) (since all clusters simulated have tpo ^ 5.5), signalling that no unique member of the King sequence can 
be found to represent the cluster. In our simulations, clusters always went out of bounds at v — » —0.6, so below we refer 
to this condition as 'u — + —0.6'. The second criterion for an immediate halt, called 'M20', is realized if a single event of 
j A M\ > 0.20M occurs which would greatly distort the cluster and invalidate our approximation of small AM/M. There 
are four other criteria which are recorded the first time they occur, but do not stop the simulation. These are: 'M10', a single 
occurence of 10 per cent mass loss; l \M\ decrease of the mass of the cluster to half its original value; '|J3', change of the 
energy of the cluster up or down by half of its original value; and '-73/M', change of the quantity E/M up or down by half 
of its original value. Survival of the cluster will be called criterion V. Note that the criteria M10, \E, §M, and ^E/M can 
happen repeatedly before v — > —0.6, M20, or s, but are only recorded the first time they occur. The two conditions \E and 
v — > —0.6 correspond most closely to the criteria used by previous investigators, but since our simulations include mass loss 
the correspondence is not exact. It is possible for a cluster to be destroyed in less than 10 10 years by M20 or v — > —0.6 even 
though iM, \E, M10, or \EjM might not have had a chance to occur yet. Moreover, M20 may occur before v — » —0.6. 
When analyzing the results, the competition among the various criteria must be kept in mind. 

This simulation ignores a number of possibly important effects. The most restrictive approximation is the use of the King 
sequence to model the evolution of the cluster. Relatively small changes in mass and energy may lead to ipo — > 0, so a normal 
cluster could have a lifetime larger than is found here. For example, a M ~ 10 6 Mq cluster could lose 30 per cent of its binding 
energy and mass, but no longer be fit well by a King model; yet you would still have a M ~ 7 x 1O 5 M0 cluster. Our treatment 
assumes that clusters become unstable to rapid dissolution when ipo — ► 0. 

To narrow the focus to the effects of halo black holes, we have neglected the disk and bulge components of the galaxy; 
destructive effects such as disk shocking and collisions with molecular clouds have also been suppressed entirely. Internal 
evolution of the cluster and non-spherical galactic tidal fields are not treated. Clusters are kept at single R g , so time dependence 
of the density of halo black holes and relative velocites along cluster orbits is neglected. Lastly, we will make fractional errors 
of order AM/M in our method. 

Nevertheless, our model advances previous attempts to constrain properties of a hypothetical population of halo black 
holes via their effects on clusters. We include mass loss which, as was shown in §^1], contributes significantly to cluster heating, 
sometimes slightly more than 'Spitzer' heating. Moreover, since the mean mass loss is comparable to the mean energy input, 
the evolution of v oc E/M 5 ^ 3 is driven by both AE and AM. Here no simplifications of the energy input are employed except 
the impulse approximation and the straight line orbit approximation. Hence, the Monte Carlo method for finding AE(b) 
includes the important 'shocking' effect of the v ■ Sv. Lastly, the correct expression for the rate of collisions in eq. (|l^) is used. 
We also note that many criteria for the cluster to be disrupted have been used in the past; the most popular is AE(t)/\E\ = 1. 
This choice is usually implemented with no regard to mass loss and evolution. Here, we test a variety of criteria for cluster 
destruction in order to find the most restrictive. 

One last technical detail which must be discussed is the value of fc max , the maximum impact parameter for our scattering 
experiment. Since the number of collisions that must be simulated oc t^, it is essential to choose as small a b max as possible 
without losing accuracy. First examine the small Mbh case in which the cluster cannot be destroyed in a single pass. In this 
slow heating limit, AE(b) oc b~ 4 outside the core and J Q max 2nbdbAE{b) converges as fc max . Nearly all of the heating results 
from penetrating encounters and fo max — rt will be a good approximation. However, in the large Mbh case there is a critical 
impact parameter b\ oc Mbh inside of which the cluster is destroyed, and we must choose fe max > fed oc A/^ 2 . An estimate for 
Mhigh, the value of Mbh above which & max must increase oc Mbii, is discussed in §R. 



8.4 Description of the Simulations 

The time necessary to realize a cluster with JV = 1000 stars and compute the energy input was ~ lsec on a Sun workstation. 
The computation time needed to map out the history for one cluster is N ev x lsec ~ (80hrs)(M/Mbh)(10pc/rt). This severely 



N' 



N-N' 




(74) 
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limited practical choices for iV and m. In order for the simulations to be realistic, we chose to simulate the full range of Mbh 
for AM 4 only, using m = 0.7M Q and N = 1000 initially. 

To study much more massive clusters, the time limitation would force the number of stars used to be a small fraction 
of the physical value. As a consequence, the variance in energy input and mass loss would be unrealistically large in the 
simulations, and hence the Green's function would be spread over too broad a range of final states. Too small a number of 
stars will decrease the number of stars in the portion of phase space from which particles are ejected resulting in an incorrect 
evaluation of mass loss, especially if the total mass lost in the simulation is small. Uneven sampling of phase space, arranged 
with finer sampling near the escape surface, could alleviate this problem. 

For a given Mbh, &max, and /halo, a number Nt trials were performed. The number of clusters destroyed, iVdest, by each 
of the six criteria was recorded and the fraction /dcst = N^est/N of clusters destroyed computed. In addition to /dest for 
each criterion, the sum of /dest(^ — * —0.6) + /dost (-^20) is computed because the sum of the fraction destroyed by these two 
conditions is not subject to competition effects which appear for the six criterion separately. 

At least two runs with different 6 ma x were done for each cluster at a certain Mbh. The runs with the larger 6 max are 
presented here. The variation in /dest for the two runs was in all cases within the error bars shown in the figures. The value of 
bmax used for the results presented here was 6max = 2r t for Afbh < 100M and b max = 2r t (M bb /W0M Q ) 1/2 for M hh > 100M. 
These values of 6 ma x are much larger than is needed, as M^igh ^4x 10 s Mq = 570M (see Table 3), and the destructive radius 
is not outside r t until Afbh > Afhigh. 

The number, Nt, of trials ranged between Nt = 40 for small Mbh and Nt = 1000 for large Mbh- The number of trials was 
restricted by the run time, which was greatly increased for small Mbh because of the large number of collisions oc /haio/Mbh- 
Bayesian methods were used to compute P s and its uncertainty. The error bars in the figures span the range of P s around the 
peak of its posterior containing 68 per cent probability. 



8.5 Results and Discussion 

The first set of runs was for /halo = 1- The fractions of trial clusters destroyed by our various criteria are plotted in Fig. |^. 

The dominant destruction mechanism by far is v — > —0.6. For this mechanism, the black hole mass at which 50 per cent 
of the clusters were destroyed is 

M b h,crit(50%) ~ 600M Q =0.86M (75) 

The other destruction mechanisms only come into play once Mbh <;5x 1O 4 M0. At this point, v — > —0.6 still destroys 100 per 
cent of the clusters (as predicted in Table 4 for /halo = 1-0); even the approximately 40 per cent that die via M20 do so in 
the same encounter that pushes v — ► —0.6. If v — > —0.6 were not stopping the simulations simulations the other destruction 
criteria would be more important, and /dest would not necessarily asymptote to small values as in Fig. ^. 

We can compare the results in Fig. [)] for the criterion v — > —0.6 with the Gaussian Monte-Carlo simulations in Fig. ^[ 
There the critical black hole mass is Mbh,crit(50%) ~ 65OM0, which is different from the full Monte-Carlo simulation by <10 
per cent. Considering that the error bars in /dest are of order 10 — 15 per cent, the two numbers agree. In addition, the mass 
range for which 0.9 > /dcst > 0.1 is <5Mbh,crit ~ 250M© for both the Gaussian Monte-Carlo and full Monte-Carlo cases. The 
agreement of these two methods shows that the assumptions of (1) a Gaussian distribution of energy input and mass loss, 
and (2) smooth Fokker-Planck evolution, are accurate in the small Mbh regime. In addition, the ansatz for the variance of 
the number of stars ejected in § |5.l| gave similar results are the full Monte-Carlo treatment, as evidenced by the comparable 
widths <5Mbh, C rit- 

The second set of runs explores the large M b h limit of §|. Two runs for AM 4 at M b h = 1000M = 7 x 1O 5 M and 
Mbh = 5000M = 3.5 x 10 6 Mq were performed as discussed in the previous section. The results are shown in Fig. [l| and[tl| 
In addition, we also tried NGC 5053 since it has a fairly large mass, M = 377OOM0, and yet it restricts /halo to be small. 
Because the run time would have been prohibitively long using the correct number of stars in this cluster, we chose iV = 1000 
for NGC 5053. As discussed previously, this should not change the values of P 3 much since we will be in the limit in which 
the collisions are quite strong, and the variances are expected to be small compared to the means. A single black hole mass 
of Mbh = lOOOAf = 3.77 x 10 7 M o was used to get the results shown in Fig. [l^ 

The values of / ha io used in Fig. [l(] and |ll] are (from left to right), / ha io = 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6. The values of 
Nd computed with eq.(|6l|) for each cluster will depend on the criterion for destruction through the values of / and C chosen. 
This is why the two curves for the two criteria v — > —0.6 and AM/M = 20 per cent do not have the same values of N4 for 
each data point of a given /haio- We have converted the values of /haio to the values of Nd shown in the graphs using eq. ( |d0[ ) 
and (|l|). The values of /halo used for NGC 5053 in Fig. [[2] are (from left to right), / ha i D = 0.01, 0.055, 0.128, 0.17, 0.273. 

As the results for AM 4 and NGC 5053 were a bit different, let us discuss AM 4 first. Notice that the full Monte-Carlo 
P B curves m Figs. and agree quite well with each other. This was predicted in §[7| from the fact that Nd, the number 
of destructive encounters, is independent of Mbh- Next notice that the curves for the v — > —0.6 destruction criterion agree 
quite well with eq. (^) which was used to draw the solid lines in the figures. The predicted P s for AM/M — 20 per cent did 
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Figure 9. Fraction of initial AM 4 clusters destroyed according to the various criteria for destruction as a function of Af^h- These curves 
were made using the full Monte Carlo method. 
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Figure 10. Comparison of full Monte Carlo simulations against theory for AM 4 with = 5000Af . The circles are the points of P s 
vs. Nd for the destruction criterion v — > —0.6. The triangles are for P s vs. Nd with the destruction criterion of 20 per cent mass loss in 
a single encounter. 



not, however, agree very well with eq. (|62|). This can be attributed to the fact that v — > —0.6 generally occurs well before any 
collision with AM/M = 20 per cent. In fact, the good agreement of the numerical results with eq. ( |67| ) for v —* —0.6 is partly 
a consequence of the relative rarity of collisions with j A M\/M — 0.2 for AM 4. 

For NGC 5053, the fraction of clusters destroyed by AM/M = 20 per cent is much larger, and hence neither curve agrees 
well with either eq. (^) or (^). For Nd J$ 1.5, since few clusters were destroyed by AM/M = 20 per cent the agreement for 
v — > —0.6 was good, but competion between the two destruction criteria is apparent for Nd > 1.5. However, note that the 
survival probability against v — > —0.6 is smaller than for | A M\/M — 0.2. Since the latter occurs only in a single catastrophic 
encounter, whereas the former also involves the cumulative effect of numerous gentle collisions, it is incorrect to ascribe cluster 
destruction at large M^h entirely to the effect of the single most destructive collision. 
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Figure 12. Comparison of full Monte Carlo simulations against theory for NGC 5053 with M^h = 1000M. The circles are the points 
of P s vs. N d for the destruction criterion u — > —0.6. The triangles are for P s vs. N d with the destruction criterion of 20 per cent mass 
loss in a single encounter. 



9 CONCLUSIONS AND PROSPECTS FOR FUTURE WORK 



Our results confirm the basic picture of cluster evolution proposed by Wielen (1988) in which there are two distinct regimes 
based on the size of the black hole mass; for M^h too small to destroy the cluster in a single collision, the evolution can be 
modelled by a smooth average over many collisions while for Mbh large enough to destroy the cluster in a single collision, the 
survival probability over 10 10 years depends both on the heating from non-destructive encounters and the stochastic effect of 
the individually destructive collisions. 

Several technical improvements have been made over previous investigations for the Mbh Afhigh limit. The cluster's 
structure was allowed to change by evolving it along a King sequence; comparison of the values of Mbh.crit in Fig. 8 and the 
last column of Table 3 show that evolution accelerates the dissolution process for Moore's weakly bound clusters, leading to 
a stricter limit. The inclusion of mass loss gave energy changes comparable to the usual 5v 2 /2 heating. Finite iV effects have 
been included, and shown to be relatively unimportant in determining Mbh.crit for the Mbh.crit <C Mhigh case since 'the range 
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of final states' is quite small when the cluster is evolved over times long enough to disrupt it. The Fokker-Planck model (§ 
for the evolution of an ensemble of initial clusters agreed closely with the full Monte-Carlo simulations (^) showing that the 
the 'slow heating' approximation is indeed an accurate representation of the evolution. Indeed, the simple model of §3.2 with 
constant diffusion coefficients given by the convenient formulas in Table 2 gives quick, analytical results accurate to within a 
factor of a few (Table 3). 

The final fate of the clusters studied in this paper was always dissolution. This occured for three reasons. First, we ignored 
internal relaxation which tends to drive a cluster to core collapse. The role of internal relaxation is currently being studied 
(Murali et al. 1998). Second, when the cluster is heated slowly by many penetrating encounters, the final fate is dissolution 
independent of cluster concentration. Third, in the large Mbh limit in which all non- destructive encounters are tidal, clusters 
with ipo < 5.5 dissolve and those with ?/>o > 5.5 core collapse. All the clusters examined here had ipo < 5.5. 

Our results for the slow-heating, Mbh -C Mhigh limit are given in Table 3, Fig. ^, and Fig. ^| The strictest limit on Mbh 
comes from the full Monte-Carlo calculation for AM 4 with Mbh.crit — 600M Q . For the v out of bounds criterion, several of 
Moore's clusters are disrupted at Mbh ~ lOOOM© and seven of nine are disrupted for Afbh ~ 6OOOA/0. For the AE/\E\ = 0.5 
criterion, six clusters die at Mbh ~ 5000 and two die for Mbh ~ 1500. To summarize, in this regime it is extremely unlikely 
that all the clusters could survive unscathed for 10 10 years if Mbh is greater than a few thousand solar masses. 

For the Mbh S> Mhi g h limit, the probability of survival P a does not tend to zero as Mbh oo, but instead asymptotes 
at a nonzero value (which does not have to be small). A simple Fokker-Planck model to determine P 3 has been developed 
(§(jj) in which both the close, destructive collisions and the distant, non-destructive collisions are included. The distribution 
of energies for an initial ensemble of clusters is shown to depend only on the expected number of destructive encounters for 
the cluster as a function of time, and the resultant P a gives good agreement with the full Monte-Carlo simulations (when 
the competition of the various survival criteria is small). Inclusion of the distant, non-destructive encounters leads to smaller 
values of P a which in turn gives much tighter limits on the allowed fraction of the halo in black holes. 

The results for the Mbh ^> Mhigh limit are presented in Table 4 and Fig. ^j, |l^, O, and |l^. For a particular cluster and a 
given value of P s , one can place limits on /halo- For P s = 0.5, Table 4 gives values of /halo for Moore's clusters ranging from 
0.02 to 0.9. In this "tidal" regime, it is unlikely that /halo > 0.3 since then P s < 0.1 for most of Moore's clusters. 

The existence for 10 years of the set of nine tenuous globular clusters studied by Moore (1993) places severe restrictions 
the allowed mass and fraction of the halo mass in massive black holes. Left unanswered in our paper is the question of whether 
Moore's clusters were initially similar in size to what we see today, of if they have been 'whittled down' to their present small 
stature. Indeed, even in isolation, the least massive of these clusters, AM 4, PAL 13, and perhaps NGC 7492 might be expected 
to have evaporated in ~ 10 9 years, suggesting that they are relatively young (Murali 1998; Murali et al. 1998; Spitzer 1987). 
To answer this question fully, one must evolve a representative population of initial clusters and attempt to reproduce the 
observed population today. The work for this project has already begun (Murali et al. 1998). Significantly tighter limits may 
result from these new investigations by the introduction of other sources of cluster evolution, such as internal relaxation, 
evaporation, and disk shocking. 
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APPENDIX A: THE CENTER OF MASS VELOCITY KICK 

For a continuous cluster with spherical mass density p(r), the velocity kick to the cluster center of mass takes on a very simple 
form. For b = be x , V rc i = V rc \e z , and projected cluster position R = i? cos^e^ + Rsvn{(f))e y , eq. ( fl8| ) gives the mean velocity 
kick 

(Aw) = 2GMbh 1 / d s xp(r)t b ~ Rcos ^ e - ~ tfsin^K 



Vtei M J ' b 2 +R 2 -2bRcos((t>) 



2GM bh 1 /* 7 * f 271 (b- Rcos(<j>))e x - Rsm(<j>)e y 



Kci MJ dRRS (R)J Q d4> b 2 +R2 _ 2bRcosW 

2GMbh 1 r dR2nRm) 0Jt^l ex 

b 

e x (Al) 



Vrel M 

2GM bh M cyl (b) 



bV Tei M 

where E(R) = 2 / (rt_fl 5 dzp(VR 2 + z 2 ) and M cyl (fe) = f* dR2irRE(R) is the mass enclosed within a cylindrical distance 
b from the center of the cluster. The b dependent function M cy \(b) /bM is zero at the cluster center, reaches a maximum inside 
the cluster and is equal to 1/6 outside the cluster since M cy \(b > r t ) = M, recovering the usual velocity kick for a structureless 
particle. 
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